;; Non-robust 2D Delaunay triangulation using the Bowyer-Watson algorithm.
;; Copyright (C) 2022  Rocks Mazama
;;
;; This program is free software: you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation, either version 3 of the License, or
;; (at your option) any later version.
;;
;; This program is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;; GNU General Public License for more details.
;;
;; You should have received a copy of the GNU General Public License
;; along with this program.  If not, see <https://www.gnu.org/licenses/>.

(ns mazama.compgeom.delaunay.bowyer-watson
  (:require [mazama.compgeom.core :refer [incircle
                                          on-or-left-of-line-segment?
                                          point-orient2d]]))

(defn- make-ccw
  "Return a sequence of the three non-colinear points in ccw order."
  [a b c]
  (case (point-orient2d a b c)
    -1 (vector a c b)
    0 nil
    1 (vector a b c)))

(defn- edges
  "Return the edges of a `tri`."
  [[a b c]]
  [[a b] [b c] [c a]])

(defn- twin
  "Return the twin of the given edge."
  [[u v]]
  [v u])

(defn- intriangle?
  "Return true iff the point `d` lies inside the closed triangle defined by `a`, `b`, and `c`.
  NOTE: The points a, b, and c must be in counterclockwise order,
  otherwise the result will be inverted.
  EXTENSION:
  One of the points a, b, or c can be nil which represents a ghost vertex at infinity."
  [[a b c] d]
  (cond
    (nil? a) (on-or-left-of-line-segment? b c d)
    (nil? b) (on-or-left-of-line-segment? c a d)
    (nil? c) (on-or-left-of-line-segment? a b d)
    :else (and (>= (point-orient2d a b d) 0)
               (>= (point-orient2d b c d) 0)
               (>= (point-orient2d c a d) 0))))

(defn- add-triangle
  "Add the triangle `tri` to the `mesh`. `tri` must be ccw."
  [mesh [a b c :as tri]]
  (assoc mesh
         [a b] tri
         [b c] tri
         [c a] tri))

(defn- add-triangles
  "Add the sequence of `triangles` to the `mesh`."
  [mesh triangles]
  (reduce add-triangle mesh triangles))

(defn- delete-triangle
  "Delete the `triangle` from the `mesh`."
  [mesh triangle]
  (apply dissoc mesh (edges triangle)))

(defn- delete-triangles
  "Delete the sequence of `triangles` from the `mesh`."
  [mesh triangles]
  (reduce delete-triangle mesh triangles))

(defn- child-node
  "Return the child-node corresponding to the `parent-edge`."
  [mesh p parent-edge]
  (let [child-edge (twin parent-edge)
        child-tri (mesh child-edge)]
    (cond
      (nil? child-tri) {:cavity-edge parent-edge}
      (>= (incircle child-tri p) 0) {:tri child-tri :excluded-edge child-edge}
      :else {:cavity-edge parent-edge})))

(defn- children-fn
  "Return the children of the `parent-node`."
  [mesh p parent-node]
  (let [parent-edges (-> (:tri parent-node)
                         edges
                         set
                         (disj (:excluded-edge parent-node)))]
    (map (partial child-node mesh p) parent-edges)))

(defn- cavity-seq
  "Return a seq of nodes that make up the cavity starting with `tri` containing point `p`.
  A node is a map with one of the following keys:
  :tri - a triangle making up part of the cavity.
  :cavity-edge - an edge that forms part of the cavity border.
  A :tri node can also have an optional property:
  :excluded-edge - an edge of the triangle not to be considered when looking for children."
  [mesh p tri]
  (tree-seq :tri (partial children-fn mesh p) {:tri tri}))

(defn- locate-point
  "Find the triangle in the `mesh` containing the point `p`."
  [mesh p]
  (let [mru (:mru mesh)
        tri0 (mesh [nil mru])
        tri1 (mesh [mru nil])]
    (cond
      (intriangle? tri0 p) tri0
      (intriangle? tri1 p) tri1)))

(defn- insert-point
  "Insert the point `p` into the `mesh`."
  [mesh p]
  (let [tri (locate-point mesh p)
        cav-seq (cavity-seq mesh p tri)
        old-cavity-triangles (->> cav-seq (filter :tri) (map :tri))
        new-cavity-triangles (->> cav-seq
                                  (filter :cavity-edge)
                                  (map :cavity-edge)
                                  (map (fn [[a b]] [a b p])))]
    (-> mesh
        (delete-triangles old-cavity-triangles)
        (add-triangles new-cavity-triangles)
        (assoc :mru p))))

(defn- insert-points
  "Insert the points in `coll` into the `mesh`."
  [mesh coll]
  (reduce insert-point mesh coll))

(defn- init
  "Return the initial mesh and any remaining points.
  The initial mesh is formed from the first n >= 2 colinear points and the first non-colinear
  point, together with the ghost triangles around the boundary."
  [S]
  (loop [S (-> S set sort) ; Remove duplicate points since they cause problems above.
         colinear []]
    (when-let [p (first S)]
      (if (< (count colinear) 2)
        (recur (next S) (conj colinear p))
        (let [[q r] colinear]
          (if (zero? (point-orient2d p q r))
            (recur (next S) (conj colinear p))
            (let [tri (make-ccw p q r)
                  [p q r] tri]
              (-> {:mru p}
                  (add-triangles [tri [r q nil] [p r nil] [q p nil]])
                  (insert-points (drop 2 colinear))
                  (assoc :mru p)        ;restore the :mru to its rightful location
                  (vector (next S))))))))))

(defn- clean
  "Clean up implementation cruft."
  [mesh]
  (dissoc mesh :mru))

(defn delaunay-begin
  "Return the initial mesh and a sequence of any remaining points.
  The initial mesh is formed from the first n >= 2 colinear points and the first non-colinear
  point, together with the ghost triangles around the boundary.
  The mesh may contain keys that are keywords. These are implementation details and are subject
  to change."
  [S]
  (init S))

(defn delaunay-increment
  "Insert the point `p` into the `mesh` and return the new mesh."
  [mesh p]
  (insert-point mesh p))

(defn delaunay-end
  "End the algorithm and return the final mesh."
  [mesh]
  (clean mesh))

(defn delaunay
  "Compute the Delaunay triangulation of `S` using the Bowyer-Watson algorithm.
  The triangulation is a map with the keys being half-edges and the values being the associated
  triangle. Each triangle is a sequence of points in counterclockwise order. The returned
  triangulation includes the ghost vertices and triangles around the boundary of the points.
  Section 3 http://people.eecs.berkeley.edu/~jrs/meshpapers/delnotes.pdf"
  [S]
  (when-let [[mesh S] (delaunay-begin S)]
    (->> S
         (reduce delaunay-increment mesh)
         (delaunay-end))))

(defn ghostbuster
  "Filter out the ghost vertices and triangles from the `mesh`."
  [mesh]
  (into {} (filter (fn [[[a b] [p q r]]] (and a b p q r)) mesh)))

(defn triangles
  "Return a sequence of triangles from the `mesh`."
  [mesh]
  (-> mesh vals set))
